home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / calc / mersenne.cal < prev    next >
Text File  |  1995-07-17  |  833b  |  45 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Perform a primality test of 2^p-1, for prime p>1.
  7.  */
  8.  
  9. define mersenne(p)
  10. {
  11.     local u, i, p_mask;
  12.  
  13.     /* firewall */
  14.     if (! isint(p))
  15.         quit "p is not an integer";
  16.  
  17.     /* two is a special case */
  18.     if (p == 2)
  19.         return 1;
  20.  
  21.     /* if p is not prime, then 2^p-1 is not prime */
  22.     if (! ptest(p,10))
  23.         return 0;
  24.  
  25.     /* calculate 2^p-1 for later mods */
  26.     p_mask = 2^p - 1;
  27.  
  28.     /* lltest: u(i+1) = u(i)^2 - 2 mod 2^p-1 */
  29.     u = 4;
  30.     for (i = 2; i < p; ++i) {
  31.         u = u^2 - 2;
  32.         u = u&p_mask + u>>p;
  33.         if (u > p_mask)
  34.             u = u&p_mask + 1;
  35.     }
  36.  
  37.     /* 2^p-1 is prime iff u(p) = 0 mod 2^p-1 */
  38.     return (u == p_mask);
  39. }
  40.  
  41. global lib_debug;
  42. if (lib_debug >= 0) {
  43.     print "mersenne(p) defined";
  44. }
  45.